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A COMPARISON OF COWELL VS POWER SERIES NUMERICAL 
INTEGRATION AS APPLIED TO ORBITAL CALCULATIONS 

C. E. Velez 
R. V. Borchers 

Program Systems Branch 
Mission & Trajectory Analysis Division 


ABSTRACT 

This paper presents the results of a comparison ci the classi- 
cal Cowell numerical integration method with the power series inte- 
gration method as applied to artificial satellite orbital calculations. 
These integration methods are described in detail, the relative ad- 
vantages and disadvantages of each are discussed, and numerical 
examples are given. The conclusion derived from this study is that 
if the orbital motion is governed by a simple, smooth force model, 
the power series method is most effective; however, for applications 
involving complex or non- smooth models, such as those occurring 
in definitive near-' Earth trajectory calculations, the Cowell technique 
is more efficient. For example, the power series method is shown 
to be 10 to 15 times more expensive than the Cowell method when 
used to integrate the motion of an Earth satellite such as GEOS-B; yet 
both methods achieve equivalent accuracies. 
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A COMPARISON OP COWELL VS POWER SERIES NUMERICAL 
INTEGRATION AS APPLIED TO ORBITAL CALCULATIONS 


INTRODUCTION 

The following is a report on a comparison study between two numerical integra- 
tion methods used for the integration of v .he equations of motion of a satellite 
and the concomitant variational equations. The first is the classical "second 
sum” Cowell method based on an integrated Newtonian interpolation polynomial 
and the second method involves a recurrent power series time expansion of the 
acceleration model itself. Both these methods are well known. The Cowell 
method is used in many of our current orbit determination programs and is de- 
rived in detail in Henrici [1], The power series method, although not as widely 
used, is also well known for its applications in celestial mechanics. For ex- 
ample, Steffensem 2] applies the method to the restricted three body problem. 

A similar, more recent analysis is given by Fehlberg [3]. The purpose of this 
study is to determine the relative effectiveness of the power series technique 
as applied to the integration of the equations of motion of an artificial satellite 
disturbed by a non-spherical primary and other effects, and also to the integra- 
tion of the variational equations required to form the state transition matrix. 


COMPARISON CRITERIA 

Generally, the simplest comparison criterion for different interpolatory or 
polynomial type numerical integration procedures is the number of function 
evaluations required to achieve a demanded accuracy for a fixed integration 
interval. This criterion is valid for all Runge-Kutta, Adams, Nordsieck, Cowell 
etc. methods, in situations where the effort required to evaluate the derivative 
dominates the computational effort of the method Hself. For applications to the 
accurate computation of satellite trajectories in the vicinity of a non- 
spherical body, this is always the case. This criterion, as opposed to direct 
measurement of computation time, enjoys the property that it does not depend 
on computer anomalies such as core speed, compiler quality, programmers 
skill, etc. Unfortunately, this criterion is not valid for the power series 
method since at each integration step, the function being integrated is evaluated 
together with its derivatives which generally significantly increases the compu- 
tational effort. Furthermore, in the case of orbital motion, these derivatives 
can be computed recurrently, making the actual computational cost per step for 
the method difficult to calculate in terms of function evaluations. For these 
reasons, the basic criterion for comparison used in this report is the amount of 
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computer time required by eac h of the methods to achieve a specified orbital 
accuracy in the position and |X)sition partials after a given |>eriod of integration. 

In each case an effort is made to show that the most efficient mode of operation 
available was used for each method. 

Hence, unavoidably, our numerical comparison results will not only reflect the 
relative merits of the integration methods, but also the computer programs used. 

To try to minimize any bias, the programs used were those written by the 
develojjers of the particular techniques as applied to the problem (references!*!] 
and ‘5 I), the problem parameters (orbit, force model, etc.) and integration 
parameters (stepsize, order, tolerances, etc.) were varied to demonstrate the 
performance of each method under different conditions, and the referenc e 
solutions used to measure accuracy were obtained by running each respective 
program in a high accuracy mode. 

Other considerations which should be made when comparing integration methods 
is the relative ease in which each method can be implemented and maintained 
to solve a specific problem. For the case of orbital integration, the important 
factors are ease of implementing the basic algorithm, core storage requirements, 
procedures for starting the process, interpolation, local error estimation and 
control, and a variable stcpshe process. Furthermore, varying the physical 
model and using "non-smooth" effects (drag, solar radiation pressure, thrust, 
etc.) in the disturbing function are frequently required in programs which 
integrate orbits, so that the performance of the integration under these condi- 
tions should also be considered. 


INTEGRATION PROGRAMS AND METHODS 

The programs used for the numerical comparisons were the stand-alone Cowell 
orbit generator from the GEOSTAR system as documented in reference 1 4 land 
the stand-alone power series orbit generator from the NAP system as docu- 
mented in reference ( 5 ]. Some of the basic characteristics of these programs 
pertinent to the comparisons are: 


GEOSTAR 

(i) Second sum ordinate Adams-Cowell predict, partial correct integration 
for the equations of motion; 

(ii) Corrector-only integration for the variational equations; 
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(iii) Variable order, variable stepsize error control with independent step 
and order control for the t juations of motion and variational equations 

integrations; 

(iv) Geopotential model includes full 15 x 15 field; 

(v) Variational equations for position partials with respect to initial state, 
geopotential coefficients, drag and radiation pressure parameters; 

(vi) Third body lunar-solar gravitational effects with internally generated 
moon-sun ephemeris. 


NAP 


(i) Power series integration of equations of motion and variational 
equations. 

(ii) Variable order, variable stepsize error control with independent order 
control for the equations of motion and variational equations integrations. 

(iii) Geopotential model includes full 15 x 15 field. 

(iv) Variational equations for position partials with respect to initial state, 
geopotential coefficients, drag, radiation pressure and mascon parameters. 

(Only state partials available in current version.) 

(v) N-body gravitational field with planetary ephemeris series initialized 
using JPL ephemeris data. 

Before presenting the numerical results obtained by using these programs, the 
following relative advantages and disadvantages are immediate concerning the 
two methods. 


ADVANTAGES OF POWER SERIES INTEGRATION 

(i) Assuming the order of the series can be made arbitrarily large, the 
integration error can be made arbitrarily small within the radius of convergence 
of the series, i.e., method enjoys the accuracy advantage of high order "multi- 
step” methods. 

(ii) An integration step in the series method does not rely on a history of 
ordinate values, as do multi step methods. This implies that power series 


3 



integration enjoys the advantages of a "single-stop" method in that it is self- 
starting and ran dynamically change steps! zc during the process without any 
significant effort. 

(iii) The local error induced by the truncated power series can be readily 
estimated, allowing a simple error control algorithm. 

(iv) Interpolation within the radius of convergence can be effected by a 
simple evaluation of the series. 


DISADVANTAGES OF POWER SERIES INTEGRATION 


(i) The method is not "general purpose" in the sense that each new problem 
to be solved requires either a new program or a significantly modified existing 
one. 


(ii) The method is frequently difficult or imjx)ssible to implement since 
the function being integrated has to be expanded in a time series. Such a series 
may not exist, and even if it does, it may not be possible to derive analytically. 

(iii) A program using the power series method is frequently not amenable 
to small changes in the function being integrated, making it expensive to main- 
tain as a general tool for research and development. 

(iv) rhe method requires that the function being integrated have high order 
continuous derivatives and hence is sensitive to small "non-smooth" anomalies. 


ADVANTAGES OF COWELL INTEGRATION 

(i) High order accuracy, limited only by numerical stability, can be 
easily attained. 

(ii) The method is "general purpose" in the sense that many problems 
could be solved with essentially the same program. 

(iii) The method is simple to implement, requiring only knowledge of the 
function being integrated. Furthermore, since it is nearly independent of this 
function, changes or modifications can be readily made. 

(iv) The method is insensitive to small "non-smooth" anomalies in the 
function being integrated. 

(v) The local error can be readily estimated. 
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DISADVANTAGES OF COWELL INTEGRATION 


(i) The method is not self-starting, requiring an inde|X3ndent procedure 
to start. 

(ii) The method requires, at each step, a ’’history” of ordinate values 
requiring an independent procedure to restart, after a step change. 

(iii) Interpolation has to 1x3 done by an independent method. 

(iv) The method is sensitive to high order numerical instabilities, effec- 
tively restricting the order of the method. 

An implication one can derive from these considerations is that if the function 
to be integrated is fairly simple, analytic or very smooth, and not subject to 
frequent modifications, then power series is probably the best technique to use, 
if all the necessary derivatives can be easily computed. A non trivial applica- 
tion which satisfies these conditions is the restricted three body problem, and 
the power series method has been used quite successfully for this case; see for 
example references [6] and 3!. 

In contrast, the application of power series to the problem of the motion of an 
artificial satellite governed by a non-spherical central body together with third 
body effects, radiation pressure and atmospheric drag forces is not as clear. 
The problem generally has a very complex formulation, is subject to discontinui 
ties and the model is constantly being modified or changed as our knowledge of 
the nature of the physical forces increases or the need for more accurate model 
ing arises. 


NUMERICAL RESULTS 

To try to determine the relative effectiveness of the two methods for this prob- 
lem, the NAP and GEOSTAR orbit generator programs were used to integrate 
the equation for the force due to a nonspherical central body, i.e., 


^ GM - B U 
x = - x + — _ , 

r 3 3x 

where _ 

R e\" 

— 1 (c^cosm\ + s™sinm\)P™(sini//), 

l_n= ? m=0 _ 



(i) 
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r radius from center of earth to satellite, 


earth’s semi- major axis, 

G M gravitational constant times mass of earth, 

\ geocentric longitude (|>ositive east), 

i i> = geocentric latitude, 

Fj (sin . ) associated Legendre polynomials. 

Differences in the way these programs handle planetary ephemeris, solar radi- 
ation, and drag did not allow the inclusion of these effects in the comparisons. 
However, it is felt that the above problem is sufficiently representative < give 
useful results. 

The following remarks are made in reference to the numerical results: 

(1) All integrations were made with both the GEOSTAR and NAP programs. 
The integrations were performed with a full 1 5th order gravity model. An 
explanation of the various columns in tables 1-6 appears as follows: In the case 
of the NAP integrations the columns from left to right have the following meanings; 

P P - number of terms used in integration of equations of motion 
and variational equations respectively; (P 2 = 0 indicates no 
variational equations are integrated). 

N - Truncation error expo no nt, with allowable truncation error 
given by 0.5 x 10 N * R \see Equation (8)); 

A b - range of integration step sizes used (power series integration 
was found to be most efficient with a variable stepsize mode); 

A x - Maximum modulus of position error, in meters, from a refer- 
ence solution; 

Ax- Maximum modulus of velocity error, in meters per second; 

AX(t) - number of digits in position partial — X - -. ) which agree with 
the reference solution; x o 

CPU- Central processing unit time in minutes. 

Similarly, in the case of the GEOSTAR integrations the first two columns from 
left to right are given by 


G 






















Fictitious Satellite — 90 Day Arc — Full 15 x 15 Potential — 360 
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45.0 Degrees M = 179.998287 Degrees 



Table 3 
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Table 5 
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Pj P 2 - represent orders in predictor corrector integration formulas 
for equations of motion and variational equations respectively, 

( P 2 “0 indicates no variational equations arc integrated) 

h - integration step in seconds. 

The remaining columns have meanings as already denoted above. 

(2) The three satellites ATS-1, GEOS-B, IMP-IV, plus a fictitious satellite were 
selected for the integration with 2, 4, and 90 day arcs. 

(3) The computations were performed on the IBM 300 model 75 and were run in a 
"sequential processing" mode so that the CPU times per run are representative 
of the calculation time required by each process. 

(4) The term "elliptic motion" indicates that the integration was performed with 
a point mass 2-lxxiy acceleration only, i.e., U in Equation (1) was set to zero. 

An immediate observation concerning these numerical results is that for the 
integration of a near-earth circular satellite such as GEOS-B (Table 4), the 
power series method was 10-15 times more expensive to use, and that for 
satellites which are further away, such as ATS-A or the fictitious satellite 
given in table 2, this ratio of power series time to Cowell time drops to 5 or 6. 
Such results indicate the sensitivity of the power series to the non-smooth per- 
turbative effects of the geopotential. This is again exemplified in tables 5 and 
6, where for the highly elliptic IMP-F case, the two body runs show power ser- 
ies superior to Cowell, but the addition of the perturbative potential increases 
the power series /Cowell ratio from less than 1 to 2 or 3. To see why the power 
series method is slower although it allows larger stepsizes, the following sec- 
tions present an analysis of the two methods as applied to the two body problem. 
Although this analysis does not reflect the actual computations performed using 
the geopotential, it details how these methods are constructed and the work re- 
quired to perform an integration step. 


POWER SERIES SOLUTION FOR TWO BODY EQUATIONS OF MOTION 
We wish to solve the two body problem 


a? G M _ 

X = - X 

1-3 

with initial conditions x Q and x () given at time t . 

The position of the satellite requires the following series with K terms 


( 2 ) 
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r (t) 



R, 0 - t 0 ) i 


K 

x (t> 22 (t ■ i o)‘ 

j “0 


( 3 ) 


where 


x 


0 



X 


0 


R 


i 





and series required for the determination of velocity are given by 


r (U rY'' j R | (t - t 0 ) 

j* 1 


-1 


K 

= 2 : 


* O) / j R (t - t„) 

I 1 

j“ I 


j-1 


(4) 


To determine the coefficients R ( , two auxiliary functions are introduced: 


e (o = — r e. (t - t 0 )i 


r J (t) 


L-j 

i*o 
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W (t) 


r (t) 
r (t) 


I 

)■ o 


7 w , -‘o) 1 


and using these functions, the coefficients for the required series are determined 
recursively as follows: 





- R , v.) • 


w 


K" 1 



W. 


K-l-j 


K- 1 



)“0 



K 

1 

(K + 2) (K i 1> 

i=o 



At each integration step there are 9K + 6 multiplications, 8K additions, and 6 
divisions to obtain the coefficients R K , R K in the power series. 

The number of operations involved in obtaining x(t) from its series representa- 
tion as given above is 3K multiplications and 3K additions. The velocity vector 
x(t) requires 3(2K - 1) multiplications and 3 (K - 1) additions. 

Hence, approximately 32 K + 6 operations are required for integrating the two 
bod’ equations of motion from time t to time t + h. The selection of K = 16 then 
gives a total of 518 operations. 
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POWER SERIES SOLUTION FOR THE VARIATIONAL 
EQUATIONS OF TWO BODY MOTION 

The acceleration for unjKjrturbed motion is given by the vector equation 


r 3 

The state transition matrix relating variations in position at time t to position 
at time t 0 is given bv the first and second fundamental matrix solutions of the 
variational system 


y = (' l + x X 1 ) ( t ) y = A ( t ) y 


(5) 


where 


- G M r ~ \ x = (|) x , <p 3 GM r ~ s , 
13 x 3 identity matrix. 

The G x G state transition is 



X(t. t 0 ) 


U 




where U (t. t 0 )and V (t, t () ) are 3 x 3 matrices satisfying (5) with initial 
conditions 


U(t 0 , t 0 ) = I. V(t # . t 0 ) = 0 
U(t 0 .t 0 ) = 0. V(t 0 , t„) = I. 


1G 


In addition to the auxiliary functions (t) and w(t), which wore introduced 
previously for the equations of motion, another auxiliary function ^ (t) given by 


9 (t) 


3 G M 

r s (t) 





J 


is used for the variational equations. 
The required series are given by 


x (t) - 9 x (t) 



K 

uo, < 0 ) = r u , (*-*o)* 

j»0 


K 

V(t. *o> *J2 V j (t " *0)* 
| a 0 


( 6 ) 


which are needed together with the series for 'x(t) given previously. 


The coefficients in these series are obtained in a recursive manner as follows 


K- I 


‘‘k = /.*i Vi-i (series for <jp) 


i *o 
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A 

■ - £ ♦, v, 

i m o 


(scries for x column matrix 


a. h k i y * r t 


K • / ( C s< * r i es tor A (i ( symme trie ma t r i x) 

i 


Uk,j= (Kl 1UK ,2) 2^ A < U «~* (SPripS f ° r U 3-3> 


K*2 


i»0 


(K 


A 

^ A V„ , (series for V, ) 

♦ 1) (K ♦ 2) Z— , J K " J 3 * 3 


i * o 


Hence we have 


X(t. t 0 ) 



1 

K 

22 u . 

j «*0 

/^v 

0 

1 

“Us]*~ 

3x3 

3 x 3 

K 

22 j u . ( » -‘ov* 1 

j 3 1 

K 

22 i v , (t-to)*' 1 

J * 1 

3x3 

3 < 3 


( 7 ) 


Each integration step is found to require approximately (>7K + G8 multiplications, 
49K + 4 additions, and 2 divisions or approximately 11GK + 74 operations 
for the U , V coefficients in the above series. 

j j 
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The number of operations required to obtain the matrices U 3i( , and V,., is 18K 
multiplications and IsK additions. The derivatives U and V require a total of 
30K - 18 multiplications and 36K - 18 additions. The sum of all the above 
operations is approximately 22 IK + 38 operations. The selection of K= 12 
yields approximately 272<’> operations for obtaining the position and velocity 
partials by the power series method. 


Power Series Variable Stepsi/.e Integration 
Consider the power series for position 


K 

x m s o - t 0 )i. 
1 *0 


Assume the modulus of the last term in the series is required to lx? equal to 
some truncation parameter * . Then the appropriate step size may be computed 
by solving 




€ 


2 


i.e., 



(«) 


In NAP, the truncation parameter e is chosen to be 0.5 x 10" N x R f where 
N is a given truncation error exponent, 

R is the modulus of x(t). 

COWELL INTEGRATION OF THE EQUATIONS OF MOTION 
Predictor- Corrector Integration Formulas 

The integration formulas used for integrating the equations of motion are given 
by [4]: 
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n* 1 


h J 

r 


I I 


s 


7 



(Stormer 


predictor) 


l, 2 


n* 1 c 


n s 


r 

L 

i *o 


a x 
i n*l“ i 


(Cowell corrector) 


x . , = h 

n* 1 p 


y 


A s , 

» n- 1 


i *0 


(Adnms predictor) 


( 9 ) 


x . , h 

n* 1 * 


■s 

n Z_ 

i*0 


fl* X 

t iv* 1“ i 


(Moulton corrector) 


where 


h the stepsize used for the integration of the equations ^f motion, 

e 

K ( = p - 2, where p is the order of the formulas, 

*S , n S the first and second sums of the accelerations, 

n n 

a., a* g ■" = the integration coefficients for the ordinate form, 

'S = V* 1 5f , 

n n 

n S = V" 1 r S = V -2 X . 

n n n 


Algorithm for Integrating the Equations of Motion 
The procedure for integrating the equations of motion is as follows: 
(i) Obtain a set of "starting” values tor the accelerations 

x K , i = 1. 2, • • • , K 

e - i 

and sums 1 S R , 1 1 S K using an independent procedure; 
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(ii) Obtain predicted values of and x' 1 ’* using the Stormer and 
Adams predictor formulas with n = K ; 

(iii) Compute the acceleration using the last predicted position and 

velocity vector; 

(iv) Obtain corrected values of and x (c ] using the Cowell and 

' ' n^l n + 1 

Moulton corrector formulas; 


(v) Compute the magnitude of the vector 



and compare with a predictor-corrector tolerance. A predictor- 
corrector cycle is completed when this vector magnitude is sufficiently 
small. A maximum number of three iterations is allowed. 


(vi) Compute new first and second sums by 




1 1 


n + 1 


i 


n + 1 



( 10 ) 


This completes one integration step and the integration is continued 
by repeating steps (ii) through (vi) with n replaced by n + 1. 


In the case in which the accelerations acting on the satellite do not depend on 
the velocity, the velocity vector is not predicted, and it is computed only after 
convergence of the position is obtained by the Cowell Corrector. 


The number of mathematical operations required in a single step of Cowell inte- 
gration of the equations of motion i "mined from Equations (9) and (11). 
Each Cowell integration step requires approximately 24 K p + 26 operations. The 
selection of K p = 9 corresponding to order 11 gives a total of 242 operations. 


The ratio for the number of operations of power series to Cowell is approxi- 
mately 518/242, or about 2/1. The introduction of perturbative effects in the 
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equations of motion will, of course, greatly increase the relative number of op- 
erations in the power series method because of the additional higher order 
derivatives which are required. These results explain, partially, the relative 
inefficiency of the power series method even though it generally allows larger 
stepsizes. 

Variable Stepsize Integration 


Stepsize modification is designed specifically for the integration of the equations 
of motion and the concomitant variational equations associated with highly ec- 
centric orbits and where fixed stepsize integration is extremely inefficient. 

The variation in stepsize is based on the concept of local error control by 
means of the following constraint equation 


where T, and T~ are specified error bounds. If the local error satisfies R > T , 

1 2 r n 1 

the stepsize is decreased, and if R n > T 2 , the stepsize is increased. The local 
error is estimated by the formula 


R ^ c 



n 



where c is a constant and x p , x c are the predicted and corrected values of the 

n 1 n 1 

position vectors. The stepsize is computed from the equation 


h 

new 




( 12 ) 


where T^ is a specified value for the local error in the range T ( < T^ < T , and 
where K is the number of backpoints used in the integration. 


COWELL INTEGRATION OF THE VARIATIONAL EQUATIONS 

The position and velocity partials as given in reference [4] for velocity free 
accelerations can be expressed as 
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(letting [U(t„. t 0 ). V(t„. t 0 )] 3 . # , [U. V]„). 


[U. v ]„ f , = [I -H]-> (W n ] 3 „ 


( 13 ) 


where 


[U.V] ntI ♦ (Y „] 3 . 6 


H = h 2 a* A , 

v 0 n + 1 


w = h 2 

n v 


I II 


V" 

L-, 

J « 1 


a* [li. V] . 

J n + 1 - j 


Y = h 

n v 


v 

L 

j = 1 


# [U. V ] 


n + 1 - j 


h v = the step size used for integrating the variational equations, 

K = p - 2, where p is the order of the integration formulas 
used for the variational equations, 

! P , ir P = che first and second sums 3x6 matrices of the acceleration 

n n 

partials, 

a*, fi* = the ordinate form of the integration coefficients corresponding 
to K . 


The inversion of the matrix (I - H) when performed by Gaussian elimination 
requires n 3 /3 + n 2 - n/3 multiplications and divisions, where n is the order of 
the matrix [7], For the case considered here, n - 3, and therefore about 17 
operations are required for the inversion. The calculation of position and 
velocity partials by Cowell integration requires approximately 72 K v + 342 
operations. For K y =5, corresponding to order 7, this gives a total of 702 
operations. The ratio for operations of power series to Cowell for the variational 
equations is approximately 2726/702 or about 4/1. 


23 



SUMMARY AND CONCLUSIONS 


The power series and Cowell methods of integration have been compared, and it 
lias been shown that although the power series is a highly stable, accurate tech- 
nique for orbital integration, it is considerably more expensive to use and main- 
tain compared to the Cowell method when applied to the problem of integrating 
highly perturbed orbital motion. On the other hand, power series integration 
affords a very accurate analytical tool which may offer significant advantages 
for error analysis work or the integration of low perturbed orbital motion such 
as transfer trajectories. 
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